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Abstract: In these lecture notes I review the theory of the non-linear evolution of 
cosmological perturbations in a self-gravitating collisionless medium, with vanishing 
vorticity. The problem is first analyzed in the context of the Newtonian approximation, 
where the basic properties of the Zel'dovich, frozen-flow and adhesion algorithms are 
introduced. An exact general relativistic formalism is then presented and it is shown 
how the Newtonian limit, both in Lagrangian and Eulerian coordinates, can be re- 
covered. A general discussion on the possible role of possible relativistic effects in the 
cosmological structure formation context is finally given. 



1 Introduction 



An important theoretical issue in cosmology is to understand the physical pro- 
cesses that occurred during the gravitational collapse of the matter which gave 
rise to the observed large-scale structure of the universe. A complementary issue 
is to reconstruct the initial conditions of the clustering process, e.g. the value 
of the cosmological parameters, the type of dark matter, the statistics of the 
primordial perturbations, starting from observational data such as the spatial 
distribution of galaxies or their peculiar velocities. Much work has been recently 
done in the latter direction, since more and more data on peculiar velocities of 
galaxies, as well as very large and complete galaxy redshift surveys have become 
available. 

A widely applied and well-motivated approximation when dealing with the 
dynamics of dark matter, either cold or hot, is to treat it as a system of particles 
having negligible non-gravitational interactions, a self-gravitating collisionless 
system. The dynamics of such a system is usually approached by different tech- 
niques, depending on the specific application. For instance, the evolution of small 
perturbations on a Friedmann-Robertson- Walker (FRW) background is followed 
by analytical methods. The non-linear evolution in cases where some symme- 
tries are present can also sometimes be followed analytically: typical examples 
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being the spherical top-hat model in the frame of the Newtonian approximation 
and the Tolman-Bondi solution in General Relativity (GR). There arc also use- 
ful approximations valid in the mildly non-linear regime, such as the Zel'dovich 
approximation (Zel'dovich 1970). Besides this classical approach, a number of 
variants have been proposed, all trying to overcome the inability to follow the 
development of structures beyond caustic formation. Different approximations 
apply in the highly non-linear regime, such as the hierarchical closure ansatz 
for the BBGKY equations (e.g. Peebles 1980) . The most general problem of 
studying the fully non-linear dynamics of a collisionlcss system in Newtonian 
theory can only be followed by numerical techniques, such as AT— body codes (e.g. 
Hockney & Eastwood 1981; Efstathiou et al. 1985). Finally, the non-linear rel- 
ativistic evolution of a pressureless fluid has been recently studied in a number 
of papers (e.g. Matarrese & Terranova 1995, and references therein). 

Here we review various methods, trying to show their possible interconnec- 
tions. Section 2 presents the Zel'dovich, frozen-flow and adhesion approxima- 
tions. Section 3 deals with the relativistic dynamics of a self-gravitating colli- 
sionless and irrotational fluid within GR. A final general discussion is given in 
Section 4. 



2 Newtonian dynamics of self-gravitating collisionless 
matter 

2.1 General formalism 

Let us start by writing the Newtonian equations for the evolution of collisionless 
particles in the expanding universe (e.g. Peebles 1980). These can be written 
using suitably rescaled variables and in comoving coordinates. We shall assume 
that the universe is spatially flat and matter dominated, so that the scale factor 
reads a(t) — ao(t/to) 2 ^ 3 (a subscript will be used to define quantities at some 
"initial time" to). A generalization of these formulae to the open (and closed) 
universe case is given in (Catelan et al. 1995). The Euler equations read 

^ U + ^ u ^ V (1) 

da 2a 2a ' 

where u = dx/da is a rescaled comoving peculiar velocity field. The symbol ^ 
stands for the total (convective) derivative 

^ ^ + V (2) 
da da 

The continuity equation can be written in terms of the comoving matter 
density i](x,t) = g(x,t) a 3 (t)/g aQ (where g a is the mean mass density at to) 



^ + ,V.u = 0, 
da 



(3) 
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while the rescaled local (or peculiar) gravitational potential tp = (3iQ/2aQ)<p s (x, t) 
is determined by local density inhomogeneities <5(x, t) = r/(x, t) — 1 through the 
Poisson equation 

W=-- (4) 

a 

We can restrict the analysis to initially irrotational flow. According to Kelvin's 
circulation theorem, in the absence of dissipation, vorticity is conserved along 
each fluid trajectory; in particular, a fluid with vanishing initial vorticity will for- 
ever remain irrotational. However, for a collisionless fluid such a property breaks 
down after caustic formation: a vorticity component is created in multi-stream 
regions, simply because the local Eulerian velocity field takes contributions from 
different Lagrangian fluid elements at the same position. Outside the regions of 
orbit mixing and/or after applying a suitable low-pass filter one can define a 
velocity potential by u(x, a) = V<£(x, a). It is then easy to derive from the Euler 
equation a Bernoulli equation relating the velocity potential to the gravita- 
tional one <p, 

By integrating along the particle trajectory x(a) one finds a formal solution 
of the continuity equation, 

7](x,a) = ?7o(q)cxp|-y da'V ■ u[x(q, a'), a']|. (6) 

This can be compared with the formula obtained from mass conservation, 
ry(x, a)d 3 x = rj (q)d 3 q, where q is the initial (Lagrangian) position of the parti- 
cle which has reached the (Eulerian) position x by the time a(t). One has either 
the well-known Lagrangian form 

7 1 [d( X ,a),a}=r l0 (ci)\\dx/dq\\-\ (7) 

where is the Jacobian determinant of the transformation x — > y, or the 

Eulerian one 

r ) (x,a)=r 10 (c l )\\dq/dx\\, (8) 

which however requires inverting the trajectory to find q(x, a). These solutions 
are only valid as long as no shell-crossing (caustic) has occurred, so that there is 
a one-to-one correspondence between Lagrangian and Eulerian positions. Before 
caustic formation all these forms are identical. 



2.2 Zel'dovich approximation 

The simplest approximation is of course the linear one, which consists in neglect- 
ing the terms u- Vu in the Euler equations and V- (Su) in the continuity equation. 
The resulting expressions read Ul/at(x, a) = — V(^o( x )> <Pl/jv(x, a) = ^o( x ) and 
Vlin(x, a ) = l+oV 2 ^o( x ), having neglected the contribution of decaying modes. 
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The next step is the Zel'dovich approximation (ZEL), based on the ansatz 
of extrapolating the equation u = — Vip (i.e. <P = —ip) beyond linear theory; 
replacing this ansatz into the Euler equations gives 

which has to be solved together with Eq.(3). The resulting system can be taken 
as the definition of ZEL. In this approximation a particle initially placed in q 
moves along a straight-line with constant "speed" u determined by the value of 
the initial peculiar gravitational potential in q, namely 

x(q,r a ) =q-T a V q(/ 2o(q), (10) 

with r a = a — a . The velocity field is conserved along each particle trajectory: 

u(x, r ) = u (q) = -Vq^o(q)- (11) 

The velocity potential, defined by u(x,r ) = V x ^(x, r a ) = V q ^o(q), which 
obeys the Hamilton-Jacobi equation 

The solution of Eq.(12) is 

*(x,r„) =<Po(q) + (x-q) 2 /2r a . (13) 

The density field is usually represented in the Lagrangian form 7?(q,r ) = 
?7o(q)||l + T a Do(q)| _ , where 1 is the unit matrix and Do the deformation ten- 
sor, with components D 0yi j(q) = d 2 <Po(q)/dqidqj. The deformation tensor can 
be locally diagonalized, by going to principal axes Xi,X 2 ,X 3 , with eigenvalues 
a\, a 2 ,a 3 . We thus have 

??( ' q ' Ta ' ) (l + T a ai(q))(l + T a a 2 (q))(l + T a a 3 (q)) ' 

According to the latter expression a singularity (caustic) in Lagrangian space 
would form at every point q where at least one eigenvalue, say a\, is negative. 

Besides being, by construction, consistent with the growing mode of linear 
perturbations at early times, ZEL provides a good approximation up to the 
time of first shell crossing. The inconsistency of ZEL can be seen as follows. 
Inserting the ansatz <P = — into the Poisson equation one gets an expression 
for the density fluctuation, Sjjyn = — aV • u, which is nothing but the linear 
theory relation between peculiar velocity and density fluctuation. This point has 
been discussed by Nusser et al. (1991), who refer to this determination of the 
density, tjdyn = 1 + $dyn, as dynamical density, to distinguish it from the 
continuity density obtained from Eq.(14). It is possible to go further this way by 
replacing this expression into Eq.(3): one gets the mass density in Lagrangian 
form i] DY N(<l,T a ) = »?o(q)/(l - r 5+(q)), where <5+(q) = -(«i(q) + a 2 (q) + 
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Q!3(q)), and the Poisson equation was used to relate ipo to So and we defined the 
(scaled) initial growing mode 5+ = 5o/ao- 

It is then clear that the Zel'dovich ansatz is only exact for one-dimensional 
perturbations, where the two above expressions for the density coincide; in 
the general three-dimensional case it fails. An alternative understanding of the 
Zel'dovich approximation might be obtained in the frame of perturbation theory 
in Lagrangian coordinates. This approach has been recently reviewed by Buchert 



2.3 Frozen flow approximation 

The frozen-flow approximation (FFA) (Matarrese et al. 1992) can be defined 
as the exact solution of the linearized Euler equations, where in the r.h.s. the 
growing mode of the linear gravitational potential is assumed. Such an equation 
is solved by u(x, r ) = u (x) = — V x (^ ( x )> phis a negligible decaying mode. In 
this approximation the peculiar velocity field u(x, a) is frozen at each point to 
its initial value, that is 



which is the condition for steady flow. Such an equation can be used, together 
with the continuity equation to define FFA. Particle trajectories in FFA are 
described by the integral equation 



particles during their motion update at each infinitesimal step their velocity to 
the local value of the linear velocity field, without memory of their previous 
motion, i.e. without inertia. This would be the case of a particle moving under 
the influence of a force in a medium with very large viscosity: the damping here 
is determined by the Hubble drag while the force is the gravitational one. 

As we shall see shortly, no caustics are formed at finite time in FFA, so that 
all expressions for the density can be used interchangeably. It is nevertheless 
interesting to write the comoving mass density as given by Eq.(6), namely 



having assumed r/o(x*) « 1. The logarithm of the density in x is given by the 
integral of the linear density field over the trajectory of the particle which has 
arrived to x at time r a , starting from the Lagrangian position q. 

FFA is, by construction, consistent with linear theory, so it follows correctly 
the evolution at early times (precisely as it happens for ZEL). The assumption 
of keeping the linear approximation for the velocity potential beyond the linear 
regime is justified by the fact that this quantity is more sensitive to large wave- 
length modes than the density, it is therefore less affected by strongly non-linear 
evolution. Stream-lines are frozen to their initial shape, so multi-stream regions 



(1995). 




(16) 




(17) 
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cannot form, unless they were already present in the initial velocity field. FFA 
therefore avoids the formation of caustics at finite r a , so one can try to extrap- 
olate the approximation after the time at which the first shell-crossing would 
have appeared according to ZEL. A particle moving according to FFA has zero 
component of the velocity in a place where the same component of the initial 
gravitational force is zero, it will then slow down its motion in that direction 
approaching such a position: particles in FFA need infinite time to reach those 
places where a pancake, a filament or a knot will occur. Moreover, since, unlike 
ZEL, these particles move along curved paths, once they come close to pancake 
configurations they curve their trajectories, moving almost parallel to them, try- 
ing to reach the position of filaments. Again they cannot cross it, so they modify 
their motion, asymptotically approaching them, to finally fall, for r a — ► oo, into 
the knots corresponding to the minima of the initial gravitational potential. Alto- 
gether, this type of dynamics implies an artificial thickening of particles around 
pancakes, filaments and knots, which mimics the real gravitational clustering 
around this type of structures. The physical thickening of the particle density 
around pancakes, filaments and knots, caused by the damped oscillations around 
these structures is replaced by an approximately exponential slowing down of 
particle motions, which however overestimates the actual particle deceleration. 
In the specular process of evacuation of initially underdense regions, FFA over- 
shoots the actual dynamics. Provided one gives up resolving the trajectories of 
individual particles, these effects produce a density field which looks roughly sim- 
ilar to the real one; in this sense the method should be considered intrinsically 
Eulerian. 

Matarrese et al. (1992) have considered the evolution of structures on large 
scales as described by FFA and compared it with the results of ZEL and of a N- 
body simulation, assuming a standard Cold Dark Matter model. Compared to 
the PM results, FFA recovers all the main structures in the correct places, even 
though they look thicker and the voids appear more empty and conspicuous. 
FFA leads to an excess of sub-structure, which is left on the way during the 
evolution instead of being erased by the hierarchical clustering process as in 
the true dynamics. The structures obtained by the Zel'dovich approximation, 
instead, are less prominent and more fuzzy, as the particles have diffused away 
from the caustic positions after shell-crossing. These results suggest that FFA is 
able to reproduce, at the statistical level, the clustering properties of the universe 
even on scales reached by the non-linear action of gravity. 

2.4 Adhesion approximation 

In the adhesion approximation (e.g. Gurbatov, Saichev & Shandarin 1989; Kof- 
man, Pogosyan & Shandarin 1990; Weinberg & Gunn 1990; Williams et al. 1991; 
Kofman et al. 1992) one modifies the Zel'dovich approach by adding an artificial 
viscosity term to Eq.(9), which is thus replaced by 



(18) 
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The viscosity is introduced to mimic the actual sticking of particles around 
pancakes, caused by the action of gravity even in a collisionless medium. The 
parameter v plays the role of a coefficient of kinematical viscosity, which controls 
the thickness of pancakes. 

The previous equation is the vector generalization of the well-known non 
linear diffusion or Burgers equation of strong turbulence (e.g. Burgers 1974). One 
can still define a velocity potential through u = V<£, which can be determined 
through the Hopf-Cole substitution <P = —2vhiU; the scalar field U satisfies the 
linear diffusion or Fokker-Planck equation, dU/dr a = ^V 2 J7, with the initial 
condition Uo(x) = exp[— #o(x)/2i/]. The resulting velocity potential reads 



<Z>(x,r a ) = -2vln 



(A-Kl>T a ) 3 / 2 



j rf 3 gcxp^-^-5(x,q,r a ; 



(19) 



where one defines the action 5(x, q, r a ) = <Po(q) + (x — q) 2 /2-r a , satisfying the 
Hamilton-Jacobi Eq.(12). The corresponding velocity field is easily obtained 
by differentiation; the Eulerian positions of the particles are found by direct 
integration of the integral equation 



rr a 

x(q,r ) = q+ / dr>[x(q,r^ 
Jo 



(20) 



while the density field can be obtained from Eq.(7). 

The Burgers equation is usually considered in the limit of small (but non 
vanishing) v, which corresponds to the limit of large Reynolds numbers, IZo = 
uoio/v, uq and £o being the characteristic amplitude and scale of the initial 
velocity field. The product uoio can be estimated either from the rms initial 
velocity potential smoothed on some scale R, (^q(R)) 1 ^ 2 , if this is convergent, 
or from the square root of the structural function of <Pq 

D{r) = <[<2> (x)-<2>o(x+r)] 2 > = (1/tt 2 ) / dkk 2 V v (k)W 2 (kR)[l-j (kr)}, (21) 

Jo 

evaluated at a suitable lag, e.g. r w R. Here V v is the power-spectrum of the 
initial gravitational potential and W(kR) a suitable low-pass filter. 

In the small v case the solution takes a simplified form which can be obtained 
from Eq.(ll) through a saddle-point approximation, 



<Z>(x,r a ) w -2i/ln 



X^faa) 1/2ex p(-7^;< S '( x >qa:> T a 



(22) 



where J(q) = ||1 + r a D (q)|, D is the deformation tensor and q a are the La- 
grangian points which minimize the action S at given x and T a . The Zel'dovich 
approximation is recovered in the limit v — > 0. This model has been applied to 
perform numerical simulations of the large-scale structure of the universe or to 
obtain some physical insight into the structure formation process in simplified 
cases. The model allows to obtain the skeleton of the large-scale matter distribu- 
tion by a geometrical technique based on the insertion of osculating paraboloids 
into the hypersurface <po(q)- 
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Let us finally mention that a detailed statistical comparison of the various 
approximation schemes discussed here against full numerical simulations has 
been recently performed by Sathyaprakash et al. (1995). 

3 Relativistic dynamics of a self-gravitating collisionless 
fluid 

3.1 General formalism 

We start by writing the Einstein's equation for a perfect fluid of irrotational 
dust. The formalism outlined in this section is discussed in greater detail in 
(Matarrese & Terranova 1995). With the purpose of studying gravitational in- 
stability in a FRW background, it is convenient to factor out the homogeneous 
and isotropic FRW expansion of the universe and perform a conformal rescaling 
of the metric with conformal factor a(t), the scale-factor of FRW models, and 
adopt the conformal time r, defined by dr = dt/a(t) (r, not to be confused with 
the variable r a of the previous section, is proportional to i 1 / 3 in the Einstein-de 
Sitter case). 

The line-element is then written in the form 

ds 2 = a 2 (r) [ - c 2 dr 2 + 7a/3 (q, T)dq a dq f} ] . (23) 

For later convenience let us fix the Lagrangian coordinates q a to have physical 
dimension of length and the conformal time variable r to have dimension of time. 
As a consequence the spatial metric 7 Q( g is dimensionless, as is the scale-factor 
a(r) which must be determined by solving the Friedmann equations for a perfect 
fluid of dust 

a'\ 8ttG 22 /r.A\ 
— Q b a 2 -Kc 2 , (24) 

2 

.2 



2 — + kc z = , (25) 

a \ a J 

where gb{r) is the background mean density. Here primes denote differentiation 
with respect to the conformal time r and k represents the curvature parameter of 
FRW models, which, because of our choice of dimensions, cannot be normalized 
as usual. So, for an Einstein-de Sitter universe k = 0, but for a closed (open) 
model one simply has k > (k < 0). Let us also note that the curvature 
parameter is related to a Newtonian squared time-scale kn through kn = kc 2 
(e.g. Coles & Lucchin 1995); in other words k is an intrinsically post-Newtonian 
quantity. 

By subtracting the isotropic Hubble-flow, one introduces a peculiar velocity- 
gradient tensor 

n = \l ai l,0 ■ (26) 

Thanks to the introduction of this tensor one can write the Einstein's equa- 
tions in a more cosmologically convenient form. The energy constraint, i.e. the 
time-time component of the Einstein's equations, reads 



Non-Linear Evolution of Cosmological Perturbations 



9 



i9 2 - i^iT + 4-0 - 
M a 



c 2 (K-6k) = l67:Ga 2 g b 5 



(27) 



where lZ a p(j) is the conformal Ricci curvature of the three-space, i.e. that 
corresponding to the metric 7^,3; for the background FRW solution 7™ 14 ' = 



(1 + f q 2 )- 2 5 afj , one has 7^( 7 



FRW\ 



2K8 a p. We also introduced the density 



contrast 5 = (g — Qb)/ Qb- 

The momentum constraint, i.e. the time-space components of the Einstein's 
equations, reads 



if 



fl\\ a 



•d 



(28) 



To avoid excessive proliferation of symbols, the double vertical bars are used here 
and in the following for covariant derivatives in the three-space with metric 7^/3. 

Finally, after replacing the density from the energy constraint and subtracting 
the background contribution, the evolution equation, coming from the space- 
space components of the Einstein's equations, becomes 



#y+2- 



m a -{n+2K)8° 



. (29) 



The Raychaudhuri equation for the evolution of the peculiar volume expan- 
sion scalar d becomes 



0' + -0 + + 4irGa 2 g b S = 



(30) 



The main advantage of this formalism is that there is only one dimcnsionlcss 
(tensor) variable in the equations, namely the spatial metric tensor 7a/ g, which 
is present with its partial time derivatives through ■d a j3 , and with its spatial 
gradients through the spatial Ricci curvature TZ a ^. The only remaining variable 
is the density contrast which can be written in the form 



S(q,r) = (l + <5o(q))[7(q,r)/7o(q)] 



-1/2 



(31) 



where 7 = det -f a p. A relevant advantage of having a single tensorial variable, for 
our purposes, is that there can be no extra powers of c hidden in the definition 
of different quantities. 

Our intuitive notion of Eulerian coordinates, involving a universal absolute 
time and globally flat spatial coordinates is intimately Newtonian; nevertheless it 
is possible to construct a local coordinates system which reproduces this picture 
for a suitable set of observers. Local Eulerian - FRW comoving - coordinates x A 
can be introduced, related to the Lagrangian ones q a via the Jacobian matrix 
with elements 



dx 



4 = 1,2,3, 



(32) 



where V A a (c[, t) is called deformation tensor. Each matrix element J A a labelled 
by the Eulerian index A can be thought as a three-vector, namely a triad, defined 
on the hypersurfaces of constant conformal time. They evolve according to 
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J A a ' = 0\J A 7 , (33) 

which also follows from the condition of parallel transport of the triads relative 
to q along the world-line of the corresponding fluid element D{aJ A a ) / Dt = 0. 

Our local Eulcrian coordinates are such that the spatial metric takes the 
Euclidean form Sab, i-e. 

7a/3 (q, r) = 6abJ A M, t)J b p (€1, t) . (34) 

Correspondingly the matter density can be rewritten in the suggestive form 

<?(q,r) -^(T)(l + (5o(q))[j7(q,r)/J (q)] _1 , (35) 

where J = detSF A a . Note that, contrary to the Newtonian case, it is generally 
impossible in GR to fix J$ = 1, as this would imply that the initial Lagrangian 
space is conformally flat, which is only possible if the initial perturbations vanish. 



3.2 Linear approximation in Lagrangian coordinates 

We are now ready to deal with the linearization of the equations obtained above. 
Let us then write the spatial metric tensor of the physical (i.e. perturbed) space- 
time in the form 

la/3 = 7a/3 + w at3 , (36) 

with 7 Q/ 3 the spatial metric of the background space - in our case the maximally 
symmetric FRW one, = j^p W ~ an d w a/3 a small perturbation. The only 
non-geometric quantity in our equations, namely the initial density contrast So, 
can be assumed to be much smaller than unity. 

As usual, one can take advantage of the maximal symmetry of the background 
FRW spatial sections to classify metric perturbations as scalars, vectors and 
tensors. One then writes 

W a p = XlaP + C|a/3 + ^ (£a|/3 + €/3\a) + ^ap , (37) 

with 

C\ a = K a a = 7r°g| Q = , (38) 

where a single vertical bar is used for covariant differentiation in the background 
three-space with metric J a f3- I n the above decomposition \ and ( represent 
scalar modes, vector modes and ir a p tensor modes (indices being raised by 
the contravariant background three-metric) . 

Before entering into the discussion of the equations for these perturbation 
modes, let us quote a result which will be also useful later. In the -d^ evolution 
equation and in the energy constraint the combination = 47^°^ — (1Z+ 2k) 8 a p 
and its trace appear. To first order in the metric perturbation one has 



(V 2 -2 K )7T a /3 + X |" +n X S a 



(39) 
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where V 2 (-) = (-)| 7 7 - Only the scalar mode \ an d the tensor modes contribute 
to the three-dimensional Ricci curvature. 

As well known, in linear theory scalar, vector and tensor modes are indepen- 
dent. The equation of motion for the tensor modes is obtained by linearizing the 
traceless part of the evolution equation. One has 

a' 

ir af3 " + 2-Tr af} ' - c 2 (V 2 - 2 K )ir afj = , (40) 

which is the equation for the free propagation of gravitational waves in a FRW 
background. The general solution of this equation is well-known (e.g. Weinberg 
1972) and will not be reported here. 

At the linear level, in the irrotational case, the two vector modes represent 
gauge modes which can be set to zero, £ a = 0. 

The two scalar modes are linked together through the momentum constraint, 
which leads to the relation x = Xo + «(C — Co)- The energy constraint gives 



(V 2 + 3k) 



87rGa 2 g b S , (41) 



-C +{^Ga 2 g h - «c 2 )(C- Co)- c 2 X0 
while the evolution equation gives 

C" + 2-C' = C 2 X. (42) 



An evolution equation only for the scalar mode ( can be obtained by com- 
bining together the evolution equation and the energy constraint; it reads 



(V 2 + 3k) 



C" + -C' - ^Ga 2 g b (( - Co) 



-87rGa 2 (? 6 (5o . (43) 



On the other hand, linearizing the solution of the continuity equation, gives 
5 = <5 -i(V 2 + 3 K )(C-Co) , (44) 



which replaced in the previous equation gives 

a 



5" + —S' -4TrGa 2 g b S = . (45) 



This is the well-known equation for linear density fluctuation, whose general 
solution can be found in (Peebles 1980). Once S(t) is known, one can easily 
obtain £ and \i which completely solves the linear problem. 

Eq.(43) above has been obtained in whole generality; one could have used 
instead the well-known residual gauge ambiguity of the synchronous coordinates 
to simplify its form. In fact, ( is determined up to a space-dependent scalar, 
which would neither contribute to the spatial curvature, nor to the velocity- 
gradient tensor. For instance, one could fix ( so that (V 2 + 3k)Co = — 2<5 , so 
that the ( evolution equation takes the same form as that for 5. 
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In order to better understand the physical meaning of the two scalar modes 
X and (, let us consider the simplest case of an Einstein-de Sitter background 
(k = 0), for which a(r) oc r 2 . By fixing the gauge so that V 2 £o = — 25 one 
obtains x(t) = xo an d 

C(r) = ^Xor 2 + B r- 3 , (46) 

where the amplitude B of the decaying mode is an arbitrary function of the spa- 
tial coordinates. Consistency with the Newtonian limit suggests xo = —j^fo, 
with ipo the initial peculiar gravitational potential, related to So through V 2 <^o = 
47rGa 2 l f?ob< l >o- One can then write 

C(r)--^or 2 + Bor- 3 . (47) 

This result clearly shows that, at the Newtonian level, the linearized metric 
is 7q/3 = 5 a fj + (\ a fj, while the perturbation mode \ is already post-Newtonian. 

These results also confirm the above conclusion that in the general GR case 
the initial Lagrangian spatial metric cannot be flat, i.e. Jq ^ 1, because of the 
initial "seed" post-Newtonian metric perturbation xo- 



3.3 Recovering the Newtonian approximation in the Lagrangian 
picture 

The Newtonian equations in Lagrangian form can be obtained from the full 
GR equations by an expansion in inverse powers of the speed of light; as a 
consequence of our gauge choice, however, no odd powers of c appear in the 
equations, which implies that the expansion parameter can be taken to be 1/c 2 . 

Let us then expand the spatial metric in a form analogous to that used in 
the linear perturbation analysis above. 

7«/3 = 7a/3 + e>(iV (48) 



To lowest order in our expansion, the evolution equation and the energy 
constraint imply that V°p = V a /s(j) — 0, and recalling that k = Kjy/c 2 , one 
gets IZ^ = lZ a p(j) = 0: in the Newtonian limit the spatial curvature identically 
vanishes. This important conclusion implies that 7 Q/ 3 can be transformed to Sab 
globally, i.e. that one can write = SAB-J A a J B fj, with integrable Jacobian 
matrix coefficients. In other words, at each time r there exist global Eulerian 
coordinates x A such that 

x(q,r)=q + S(q,r), (49) 

where S(q, r) is called the displacement vector, and the deformation tensor be- 
comes in this limit 
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The Newtonian Lagrangian metric can therefore be written in the form 

, \ r ( zA dS A (q,T)\f :B dS B {q,T)\ 
7a/j(q, r) = 5ab [S A a + ™ '- ) [S% + ™ >- ) . (51) 

One can rephrase the above result as follows: the Lagrangian spatial met- 
ric in the Newtonian limit is that of Euclidean three space in time dependent 
curvilinear coordinates q a , defined at each time r in terms of the Eulerian ones 
x A by inversion. As a consequence, the Christoffel symbols involved in spatial 
covariant derivatives (which will be indicated by a single bar or by a nabla op- 
erator followed by greek indices) do not vanish, but the vanishing of the spatial 
curvature implies that these covariant derivatives always commute. 

Contrary to the evolution equation and the energy constraint, the Raychaud- 
huri equation and the momentum constraint contain no explicit powers of c, and 
therefore preserve their form in going to the Newtonian limit. These equations 
therefore determine the background Newtonian metric 7^/3, i.e. they govern the 
evolution of the displacement vector S. 

The Raychaudhuri equation becomes the master equation for the Newtonian 
evolution; it takes the form 

& + + + 4nGa 2 g b (T 1/2 - 1) = , (52) 

where 

^ = \l ai % , (53) 

and, for simplicity, So = was assumed (a restriction which is, however, not at all 
mandatory). We also used the residual gauge freedom of our coordinate system 
to set 7^/3(70) = 8 a p, implying So = 1, i.e. to make Lagrangian and Eulerian 
coordinates coincide at the initial time. That this choice is indeed possible in 
the Newtonian limit can be understood from our previous linear analysis, where 
this is achieved by taking, e.g., Co = 0. 

The momentum constraint, = &,v, is actually related to the irrotation- 
ality assumption. 

Let us also notice a general property of our expression for the Lagrangian 
metric: at each time r it can be diagonalized by going to the local and instanta- 
neous principal axes of the deformation tensor. Calling %, the eigenvalues of the 
metric tensor, J a those of the Jacobian and d a those of the deformation tensor, 
one has 

7«(q,T) =J Q 2 (q,r) = (1 + J a (q,r)) 2 . (54) 

From this expression it becomes evident that, at shell-crossing, where some of 
the Jacobian eigenvalues go to zero, the related covariant metric eigenvalues just 
vanish. On the other hand, other quantities, like the matter density, the peculiar 
volume expansion scalar and some eigenvalues of the shear and tidal tensor 
will generally diverge at the location of the caustics. This diverging behaviour 
makes the description of the system extremely involved after this event. Although 
dealing with this problem is far outside the aim of the present notes, let us just 
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mention that a number of ways out are available. One can convolve the various 
dynamical variables by a suitable low-pass filter, cither at the initial time, in 
order to postpone the occurrence of shell-crossing singularities, or at the time 
when they form, in order to smooth the singular behaviour; alternatively one 
can abandon the perfect fluid picture and resort to a discrete point-like particle 
set, which automatically eliminates the possible occurrence of caustics, at least 
for generic initial data. At this level, anyway, we prefer to take a conservative 
point of view and assume that the actual range of validity of this formalism is 
up to shell-crossing. 



3.4 Recovering the Newtonian approximation in the Eulerian picture 

As demonstrated above, it is always possible, in the frame of the Newtonian 
approximation, to define a global Eulerian picture. This will be the picture of 
the fluid evolution as given by an observer that, at the point x = q + S(q, r) 
and at the time r observes the fluid moving with physical peculiar three-velocity 
v = dS/dr. From the point of view of a Lagrangian observer, who is comoving 
with the fluid, the Eulerian observer, which is located at constant x, is moving 
with three-velocity <fq(x, r)/dr = —v. 

The line-element characterizing the Newtonian approximation in the Eule- 
rian frame is well-known (e.g. Peebles 1980) 



ds 2 =a 2 (r) 



1 + 2yg( f' Tr ) c 2 dr 2 + 5 AB dx A dx l 



(55) 



with ifig the peculiar gravitational potential, determined by the mass distribution 
through the Eulerian Poisson equation, 

V 2 ^(x, r) - 4^Ga 2 (r)^(r)(5(x, r) , (56) 

where the Laplacian V 2 ,, as well as the nabla operator V, have their standard 
Euclidean meaning. The perturbation in the time-time component of the metric 
tensor here comes from the different proper time of the Eulerian and Lagrangian 
observers. 

It is now crucial to realize that all the dynamical equations obtained so far, 
being entirely expressed in terms of three-tensors, keep their form in going to 
the Eulerian picture, only provided the convective time derivatives of tensors of 
any rank (scalars, vectors and tensors) are modified as follows: 

D d dS 

This follows from the fact that, for the metric above, P^ B = P A B = P BC = 0, 
which also obviously implies that covariant derivatives with respect to x A reduce 
to partial ones. 

The irrotationality assumption now has the obvious consequence that we can 
define an Eulerian velocity potential <P V through 
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v(x,t)=V<Mx,t). (58) 
The Newtonian peculiar velocity-gradient tensor then becomes 

*» - • m 

because of which the momentum constraint gets trivially satisfied and the mag- 
netic Weyl tensor becomes identically zero in the Newtonian limit. 

We can now write the Raychaudhuri equation for the Eulerian peculiar vol- 
ume expansion scalar and use the Poisson equation to get, as a first spatial 
integral, the Euler equation 

a' 

v ' +v • Vv+ — v = -V(f g . (60) 
This can be further integrated to give the Bernoulli equation 

K + -*v + \{y$ v f = -<p g . (6i) 

All these equations would of course recover the form of Section 2, if the time 
variable and the peculiar velocity were rescaled as described at the beginning 
of that section. Having shown the equivalence of this method, in the Newtonian 
limit, with the standard one, it would be also trivial to recover the Zel'dovich 
approximation in this frame. This point is further discussed by Matarrese & 
Terranova (1995). 



4 Conclusions and discussion 



Aim of these notes was to introduce the reader to the theory of the dynamics of 
cosmological perturbations beyond the linear approximation. Let me spend this 
final section to discuss what I consider an open issue in this field. The issue is 
whether there exists a range of scales where relativistic effects and non-linear 
evolution both come into play. The standard Newtonian paradigm states that 
the lowest scale at which the approximation can be reasonably applied is set by 
the amplitude of the gravitational potential and is given by the Schwarzschild 
radius of the collapsing body, which is negligibly small for any relevant cosmo- 
logical mass scale. What is completely missing in this criterion is the role of the 
shear, which causes the presence of non-scalar contributions to the metric per- 
turbations. A non-vanishing shear component is, in fact, an unavoidable feature 
of realistic cosmological perturbations and affects the dynamics in at least three 
ways, all related to non-local effects, i.e. to the interaction of a given fluid element 
with the environment. First, at the lowest perturbative order the shear is related 
to the tidal field generated by the surrounding material by a simple proportion- 
ality law. Second, it is related to a dynamical tidal induction: the modification of 
the environment forces the fluid element to modify its shape and density. Third, 
and most important here, a non vanishing shear field leads to the generation 
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of a traceless and divergenceless metric perturbation which can be understood 
as gravitational radiation emitted by non-linear perturbations. Note that the 
two latter effects are only detected if one allows for non-scalar perturbations in 
physical quantities. Truly tensor perturbations are in fact dynamically generated 
by the gravitational instability of initially scalar perturbations, independently 
of the initial presence of gravitational waves. 

Using a post-Newtonian expansion in Lagrangian coordinates, Matarrese & 
Terranova (1995) obtained a general formula for the tensor modes ttab pro- 
duced by non-linear evolving perturbations. In the standard case, where the 
cosmological perturbations form a homogeneous and isotropic random field, 
they obtained a heuristic perturbative estimate of their amplitude in terms of 
the rms density contrast and of the ratio of the typical perturbation scale A 
to the Hubble radius th = cH -1 (where H is Hubble's constant). They found 
ftrms/c 2 ~ S 2 ms (\/rH) 2 ■ These tensor modes give rise to a stochastic background 
of gravitational waves which gets a non-negligible amplitude in the so-called 
extremely-low-frequency band (e.g. Thorne 1995), around 10 -14 - 10 -15 Hz. One 
can roughly estimate that the present-day closure density of this gravitational - 
wave background would be fi gw (\) ~ i 4 ms (V r fl) 2 ' I n standard scenarios for 
the formation of structure in the universe, the closure density on scales 1 — 10 
Mpc would be Q gw ~ 10~ 5 — 10~ 6 . The amplitude of this gravitational-wave 
contribution, ir ~ <5 2 (A/r#) 2 , is an important counter-example to the standard 
paradigm stated above, according to which relativistic effects should be propor- 
tional to ifig/c 2 ~ S(X/r H ) 2 - 
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